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Recent developments of the theoretical investigations on the one-dimensional Kondo lattice model 
by using the density matrix renormalization group (DMRG) method are discussed in this review. 
Short summaries are given for the zero-temperature DMRG, the finite-temperature DMRG, and 
also its application to dynamic quantities. 

Away from half-filling, the paramagnetic metallic state is shown to be a Tomonaga-Luttinger liquid 
with the large Fermi surface. For the large Fermi surface its size is determined by the sum of 
the densities of the conduction electrons and the localized spins. The correlation exponent K p of 
this metallic phase is smaller than 1/2. At half-filling the ground state is insulating. Excitation 
gaps are different depending on channels, the spin gap, the charge gap and the quasiparticle gap. 
Temperature dependence of the spin and charge susceptibilities and specific heat are discussed. 
Particularly interesting is the temperature dependence of various excitation spectra, which show 
unusual properties of the Kondo insulators. 



I. INTRODUCTION 

In a degenerate Fermi gas, low temperature specific 
heat is linear in T and the proportionality constant is 
given by the density of states at the Fermi energy. For 
the free electrons the density of states including the two 
spin directions is given by 



D(e F ) 



mk-p 



(1) 



where m is the free electron mass and kp is the Fermi 
wave number. 

To include the effect of electron-electron interaction, 
Landau developed the theory of Fermi liquids. Since the 
volume of the Fermi sphere is determined by the density 
of electrons, it is a reasonable assumption that the Fermi 
wave number does not change by the interaction, which 
was in fact proven later by Luttinger jlj. According to 
Landau, the effect of the interaction can be taken into 
account by replacing the bare electron mass m by an 
effective mass m*. Thus importance of the interaction 
effects in each material may be judged from the electronic 
specific heat at low temperatures. 

In ordinary metals the coefficient of the T-linear term, 
7, is of the order of mJ/mol K 2 . However, in some rare 
earth and actinide compounds there is a group of com- 
pounds whose 7 are in the range from 0.1 J/mol K 2 to 
more than 1 J/mol K 2 . This class of materials are called 
heavy Fermion systems or heavy electron systems. A 
key feature of the heavy Fermion systems is that it con- 
tains two different types of electrons: relatively localized 
/ electrons and extended conduction electrons. Interplay 
between the two degrees of freedom is an essence of the 
heavy Fermion physics. 



Since the Coulomb interaction between the / electrons 
is strong, a partially filled / shell in an isolated ion pos- 
sesses a well defined magnetic moment corresponding to 
the total angular momentum of the / shell. A weak hy- 
bridization between the / electrons and the conduction 
electrons is the source of interesting manybody problems. 

When we consider a single / shell in the sea of conduc- 
tion electrons, the magnetic moment of the / electrons is 
unstable, leading to the Kondo singlet which is a bound 
state of the / moment with the spin polarization of the 
conduction electrons pj. When we consider two / shells, 
spin polarization of the conduction electrons induced by 
an / moment tends to stabilize the magnetic moment of 
the other / shell. This is the origin of the Ruderman- 
Kittel-Kasuya-Yosida (RKKY) interaction j|. Thus the 
Kondo effect and the RKKY interaction in many cases 
compete with each other. 

If the Kondo effect dominates over the RKKY interac- 
tion by some reason, a paramagnetic heavy Fermion state 
will be stabilized. In this regime, the Kondo temperature 
or the effective Fermi temperature in a lattice problem 
sets a small energy scale at low temperatures. Existence 
of the small energy scale naturally leads to the large spe- 
cific heat since the entropy associated with the magnetic 
degrees of freedom of / orbitals should be released in the 
small temperature range. 

The most simple theoretical model for the heavy 
Fermion physics is the Kondo lattice model. The Kondo 
lattice model is given by 



H 



1 . 



(2) 



where a — (cr x ,a v ,cr z ) are the Pauli matrices and Si — 
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J2 S s> \&ss> fi S fis' is the /-electron spin at site i. In 
this review we will consider the model with only near- 
est neighbor hoppings, tij — —t for the nearest neighbor 
pairs. 

Much effort has been invested for the study of the 
model and a significant progress has been achieved in 
one dimension in the last ten years When we fix a 
lattice structure, the Kondo lattice model has only two 
parameters: one is the density of conduction electrons 
n c and the other is the strength of the exchange coupling 
normalized by the hopping energy J/t. In one dimension, 
the ground-state phase diagram in the parameter space 
is completed and shown in Fig. [I]. 

There are three different phases in the phase diagram. 
In the region from the low density limit to the strong 
coupling limit, a ferromagnetic metallic phase is stabi- 
lized. The spin quantum number in this phase is given 
by 5tot = (L — Ay/2, where L is the number of lattice 
sites and N c the number of conduction electrons. Note 
that the magnetic moment vanishes as the half-filling is 
approached. The line of half-filling is special in the sense 
that the ground state is always a nonmagnetic insulator. 
Since the lowest excitation in this phase is a spin-triplet 
excitation with a finite excitation gap, this phase is called 
an incompressible spin liquid phase. In the remaining 
part of the phase diagram (Fig.[l]) which extends from 
the weak coupling limit towards the line of half-filling, 
the ground state is metallic and paramagnetic. In this 
review we will discuss properties of the spin liquid phase 
and the paramagnetic metallic phase. 




ing 



FIG. 1. The ground-state phase diagram of the one di- 
mensional Kondo lattice model with the nearest neighbor hop- 
pings. 

Existence of the small energy scale in the Kondo lattice 
model means that correlation lengths are generally long. 
For example, to study the spin gap and the charge gap in 
the spin liquid phase the exact diagonalization was used 
at first H . However the largest system size which can be 
diagonalized by the Lanczos algorithm is just ten sites. 
This is the reason why a nontrivial form of a finite size 
scaling was necessary to obtain the functional form of the 
spin gap. 

Recently Steven White developed the density matrix 



renormalization group (DMRG) method to study the 
ground-state properties of one dimensional many body 
systems 0. Advantage of this method is that an order 
of magnitude bigger systems can be studied compared 
with the numerical exact diagonalization by the Lanczos 
algorithm. Compared with quantum Monte Carlo simu- 
lations the DMRG is free from statistical errors. The nu- 
merical errors in the DMRG come from truncation errors 
but they can be estimated from the largest eigenvalue of 
the density matrix which is truncated out. The trunca- 
tion error may be improved by increasing the number of 
basis states for the density matrix. The DMRG method is 
an ideal tool to study the one-dimensional Kondo lattice 
model and in this review we will discuss recent develop- 
ments on this subject. 

Further development of the DMRG method was 
achieved lastyear when one of the authors Q] and Wang 
and Xiang m independently succeeded to obtain ther- 
modynamic properties of the one-dimensional quantum 
XXZ model by applying the DMRG to the transfer ma- 
trix (finite-T DMRG). Application of the finite-T DMRG 
to a system with Fermion degrees of freedom started from 
the Kondo lattice model |J. 

The present article is organized as follows. In the next 
section, after a brief summary of the DMRG method for 
the ground state, we will describe the method to calcu- 
late thermodynamic properties by the finite- T DMRG 
and then extend discussions to the dynamic quantities at 
finite temperatures. In Section III nature of the para- 
magnetic metallic phase away from half-filling is shown 
to be a Tomonaga-Luttinger liquid with a large Fermi 
surface. The large Fermi surface means that the volume 
inside the Fermi surface is determined not only by the 
density of conduction electrons but also includes the lo- 
calized spins. Section IV is devoted to the discussion of 
the spin liquid phase at half-filling. After the discussion 
of the spin gap and the charge gap at zero temperature, 
we will discuss how the excitation gaps develop as the 
temperature is lowered. We will conclude the present 
review by summary and discussions in Section V. 

II. THE DENSITY MATRIX 
RENORMALIZATION GROUP METHOD 

The density matrix renormalization group (DMRG) 
method is relatively new || among the various numer- 
ical algorithms to treat many-body problems. However 
it is now widely used as one of the most standard nu- 
merical methods for low dimensional many-body sys- 
tems. In this section we first briefly outline the algorithm 
of the zero-temperature DMRG that was developed to 
study ground-state and low energy excitations of one- 
dimensional systems. The application of this method to 
the quantum transfer matrix enables us to obtain ther- 
modynamic quantities [j^]8| and the dynamical correla- 
tion functions at finite temperatures. In the second part 
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of this section we will summarize the algorithm of the 
fmitc-T DMRG. 



A. Zero temperature algorithm 

The zero-temperature DMRG method is designed to 
obtain the ground-state wave function and the low energy 
excitations with small systematic errors. The ground- 
state wave function and the low energy excitations of 
long systems are obtained by expanding the system size 
iteratively as shown in Fig. |[ The expansion of the 
system is done by putting additional sites in its central 
region to minimize undesirable boundary effects on the 
added sites. The algorithm is described in the following. 

Let us start from a system of four identical sites, for 
example, a four-sites spin chain under the open bound- 
ary conditions. An operator on the nth site, e.g. S n , is 
represented in terms of the complete basis states \i n ) as 



To increase the size of system we introduce two new sites 
between the blocks n = 1,2 and n = 3,4. Using the 
basis states |i2') and |«3'), for the new sites, we con- 
struct the Hamiltonian matrix of the expanded system 



,/3,a' 



Renaming the indices as 

ix,i 2 ' — ► 12,13' — ► k,0 ■ 



(9) 



(hi I S n 



(3) 



we repeat the procedures from the diagonalization of the 
Hamiltonian matrix. 

The key feature of the above renormalization proce- 
dure is that the new basis states \a) or |/3) of each block 
contain the information that the block is a part of the 
total system. As shown in Fig. || the edge part of each 
block connecting to the remaining part of the system is 
located in the middle of the system, and this part is not so 
sensitive to the boundary conditions imposed on the to- 
tal system. Thus we expect that the new basis states also 
dominantly contribute to the ground-state wave function 
of the expanded system which has two additional sites in 
the middle of the two blocks. 



Then we construct a representation of the Hamiltonian 
Hiiiiisii^i'^i'^ f° r the total system. The ground-state 
eigenvector 



(4) 
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is obtained by diagonalizing the Hamiltonian matrix by 
some method like the Lanczos algorithm. Then ^j 1 j 2 .; 3 .j 4 
is used to construct the density matrix 



P 



2l«2 ,l±l 2 



<3»4 



ill U/ '. , 



(5) 



for the block containing the sites n = 1 and 2. The 
density matrix specifies to what extent the basis states 
|*i)|*2) °f the block are contributing to the total wave 
function i^i^i) ■ This matrix is numerically diagonal- 
ized, and we obtain its eigenvalues A" and eigenvectors 
vf ii2 . Then we select the eigenvectors of the largest m 
eigenvalues as new basis states for the block. Here m is 
the number of the basis states kept for the block at the 
next step. Using the selected eigenvectors of the den- 
sity matrix we represent the operators on the site, for 
example, n = 2 as 



(S: 



2)a,a' 



(v a ■ )*v a - 
! \ u l 1 l 2 J 111 



(6) 



A similar procedure is repeated for the block of the sites 
71 = 3 and 4, and all operators in the original system are 
represented in terms of the new basis states 



(7) 
(8) 



L=8 



FIG. 2. Schematic diagram of the infinite system algo- 
rithm of the DMRG. 

The above algorithm is called the infinite system algo- 
rithm of the DMRG. Using this algorithm we increase the 
size of the system. In order to improve the basis states 
of the blocks, it is necessary to fix the size of system and 
use the following algorithm which is known as the finite 
system algorithm of the DMRG. The schematic diagram 
of the finite system algorithm is shown in Fig. [|. 

Let us take a block of size n — 1 on the left and an- 
other block of size n — 1 on the right whose basis states 



are represented by 



J L{r 



and 



J R(r. 



These basis 



states and representations of the operators in the blocks 
are obtained after the (n— 2)th renormalization step from 
the initial four-site system. The system of size 2n is con- 
structed by inserting additional two sites for which the 
basis states are represented by and \iw)- 

We again diagonalize the Hamiltonian matrix of this 
2n-site system H ir . ;* v ;* >* and 

obtain the ground-state 



wave function \I/ 
density matrix 



«t(n-l)»z,'»R'tR(»i 



Then construct the 



*R'ifl(»-i) 



(10) 
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for the subspace spanned by \v % L , n _ x -^j and We use 

the m important eigenvectors of this density matrix as 
a new basis states of the block containing n sites, and 
represent all operators in the new block in terms of the 
new basis states. 

In the next step, we take the new left block with 
n sites and the right block with n — 2 sites. The 
right block with n — 2 sites has been obtained at the 
(n — 3)th renormalization step from the initial four-site 
system. Inserting two sites between these blocks, we 
construct the Hamiltonian matrix of the total system 
Hi., .j ,4 . a 1 it i> i' , and repeat above all 

procedures. 

We continue to increase (decrease) the size of the left 
(right) block until the right block is reduced to the single 
site. Then we turn to decrease (increase) the size of the 
left (right) block in order to improve the basis states of 
the right block. We continue to decrease (increase) the 
size of the left (right) block until the left block is reduced 
to the single site. These procedures are continued back 
and forth until we get a good convergence. 
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FIG. 3. Schematic diagram of the finite system algorithm 
of the DMRG. 

In general, the total energy of the system is lowered as 
the basis states of the blocks are reconstructed. Thus, 
the lowest energy is obtained after the convergence. The 
wave function obtained by the finite system algorithm 
of the DMRG method may be represented by a matrix 
product form |icj . Therefore the finite system algorithm 
of the DMRG method is considered to be a numerical 
variational method which uses the matrix product wave 
function. This is another reason why we can get remark- 
able accuracy by the DMRG method. The accuracy of 
the ground-state energy and the wave function is deter- 
mined by the eigenvalues of the density matrix which 
are truncated out. Thus we can improve the accuracy 
by increasing the number of basis states m used in the 
calculations so long as the memory of computer allows. 



B. Finite temperature algorithm 

It is also possible to discuss thermodynamic quantities 
by the DMRG method, finite-T DMRG. The readers who 



are not interested in the detail of the iteration procedures 
may skip the paragraph including Eqs. ( |l3| ) to (p6|). In 
this method we use the quantum transfer matrix defined 
as 

Tn(M) = [ e -^-i,WM e -/3ft2„,2„+i/M]M 

= E E E-E E 

0-2n n <T2„tx CT2nr 2 CT 2 nT M CT2nTj£ 

(o'2n-l,TiO'2n-l,ri|e _ ^ /l2n ~ 1,2 "' M |o'2n,Ti ^2n,n) 
(o'2n,TiO'2n,T2 |e~ /3 ' l2 ™' 2 '*+ 1 / M |<72 n +l,Ti &2n+l,T 2 } 
{a2n-l,T 2 <J2n-l,T 2 \e'~ Ph2n - 1 ' 2n/M \(J2n,T 2 02n,T 2 ) 



\C2n,T M °"2n,n |C 



-0h 2 



l/M 



ff2n+l,Tjff!ln+l,n)' 

(11) 

This quantum transfer matrix is graphically shown in 
Fig. [|. Here M is the Trotter number and is the dis- 
cretized imaginary time whose intervals n+i — r,; — (3o = 
/3/M. In Eq. ( [ll|) <J2n,Ti represents states of the site 2n 
corresponding to a given imaginary time Tj . The Hamil- 
tonian H is assumed to be decomposed into two parts 

#odd = J2n=l h-2n-l,2n and H cvcn = J^t=l ^2n,2n+l SO 

that we can evaluate the matrix element of the exponen- 
tial function. Since the partition function Z is given by 
the trace of the product of the quantum transfer matrix 

Z = Tr e~^ H 

= lim Tr ( e ^ f3H ° dd/M e ~PH cvcn /M^M 



N/2 



lim Tr [TT( e -/"*»-WJ* 

M^oo - LJ - 
n=l 

' L/2 



N/2 



lim Tr 

M— >oo 



n t "( m ) 



n=l 



■Ph 2n , 2n +i/Mv\M 



(12) 



thermodynamic properties of L — > oo are determined by 
the maximum eigenvalue and its eigenvectors. To obtain 
the eigenvalue and the eigenvectors for a large Trotter 
number M, we iteratively increase the size of the quan- 
tum transfer matrix using a similar algorithm to the zero- 
temperature DMRG. 

We first represent the quantum transfer matrix as 



(M) 



n-A n-B 

(M) (M) ' 
n-A' n-B' 

(M) (M) ' 



for M : even 
for M : odd 



(13) 



Thus the transfer matrix for M — 2 is 

T{M=2) (0"2 n— l,Ti &2n— l,ri &2n— l,ra &2n— 1,T2 5 



E 



E 



n-A 

1 (M=2) 



n— l,n &2n — 1,ti C2n,T 2 j &2n,Ti &2n+\,Ti &2n+l,T 2 ) 
%M=2) ( <T 2n-l,T 2 0'2n-l,T^O'2n,-ri ! C2n,T 2 cr 2n-|-l,T2 Cr 2n+l 1 Ti)' 

(14) 
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FIG. 4. The quantum transfer matrix for M = 3. The 
^odd and ftcvcn represent /i 2 „-i,2n and /i2n,2n+i, respectively. 
Ar = Ti+i - Ti = Po and /3 = M(3 . 

Here we have introduced / ^Xf =2 ) ano - ^(M=2) which are 
denned as 

^(M=2) (°"2n— l,Ti C2n-l,ri_C r 2n,r2 I &2n,T 1 &2n+l,TiC r 2n+l,T2 ) 
= ^ (o'2n-l : TiO'2n-l,Ti|e~ /3 ' l ° £!<i ^ M |cr2 n!T1 0'2„ iTi _) 



(cr2n,Tif 2 n iT2 |e 



-/3h Bvcrl /M 



^2n+l,ri ^"2n+l,T2 

>, (15) 



-B 

'(M=2)< 

= (c r 2ri-l,T2 Cr 2n-l,T2 l e 

f 2 Tl , 7"2 



1~(M=2) { <T 2n-l,T 2 <T 2n~l,T2 Cr 2n,T 1 ^2n,T 2 <7 2n+l,T2 l7 2n+l,T 1 ) 

\ u 2n,T2°2n,T2/ 



{&2iuT2< T 2iuT 1 \e l3hcve ' l ^ M \o'2n+l,T2^27i+l,T 1 }- (16) 

where h odd = h 2n -i^n and /i OVO n = h 2n ^n+i- Then we 
iteratively increase M of ^n^n and as 



£ (M) 



2 (M) e 


^ "<(M- 


-i) 


(17) 


e 1 (M) 


^ (iWH 


-i)' 


(18) 


n-A' -/3 /i ovon 
I (M) e 


^ (AfH 


-i) 


(19) 


e 1 (M) 






(20) 



The example for the increase of M, Eqs. ( pL T|) and (18), 
for M = 2 is graphically shown in Fig. ||. 

In order to represent the transfer matrix in a restricted 
number of basis states, we have to select important ba- 
sis states which have significant weight for the represen- 
tation of the transfer matrix. For this purpose we use 
the generalized asymmetric density matrix similar to that 
used in the zero-temperature DMRG. For example, the 
density matrix which we use in the procedure Eq. (|T~ 
for M = 2 is 



p(c2n-l,TiCT2n-l,T 2 &2n+l,Ti&2n+l,T2 ) 
= E E V L ((T Tl <72n-l,T L , ^n-l^&rj 



T l " T 2 
rR 



V i(7 ri IJ2n+l,T 1 , &2n+l,T2 (7 T 2 ) 



(21) 



where V L and V R are the left and right eigenvectors of 
T(M=2) which have the maximum eigenvalue. The V L 
and V R are generally different owing to the non-Hermite 
property of the transfer matrix. The diagonalization 
of the density matrix provides eigenvectors, and v R , 
which satisfies the equations 

E V a <7t 2 )p(<7r L CrT2 \ ^'r^ ) = 7a ^ (^i <4 ) ( 22 ) 

E PWr l a T 2 ^r L (rT^v R (a TL a n )=^ ol v R {a' T y T2 ). (23) 



We select the m eigenvectors which have the largest 
eigenvalues 7 Q , and we use them as the new basis states 



\Ot2n+ll 



E E 

C"2n-l,-r 1 0"2n — l,r 2 

I, (24) 

E E 

^ Q (^2n+l,n 0"2n+l,T 2 )|°2n+l,Ti 0"2n+l 5 T- 2 

>. (25) 



Then we represent X 



A' 

(M=3) 



MS 



'^(M=3)( cr 2n-l,ria2n-l,ri,2 cr 2n-l,r2; °"2n,ri a 2n+l,Ti a^n.ra) 

= E E E~ E ~E 

&2n — 1,T1 C2n— 1,T 2 &2n+l,T 1 <72n+l,T 2 <72n,T 2 



V a (c r 2n-l,TiO'2n-l,T2 

(J 



^(M=2) (°"2n— 1,ti C2n-l,Ti cr 2ri,r 2 ; C27l,Ti C2n+l,Ti 02n+l 



{^2n-l,T 2 ^2n-l,T 2 |e P h ° dd / M \<T2n,T 2 &2iuT 2 ) 



t 2 J 

(26) 



We repeat similar procedure for Eqs. (|l7j ) to ( ^6| ) and 
obtain the maximun eigenvalue and its eigenvectros of 
the quantum transfer matrix for a desired M. 

Compared with the zero-temperature DMRG algo- 
rithm the finite-T DMRG is more subtle from the point 
of view of numerical stability. The asymmetric density 
matrix sometimes yields complex eigenvalues, although 
in principle they must be real. To avoid those unphysical 
complex eigenvalues we need accurate numerical calcula- 
tions taking account of various symmetries of the system. 

The free energy per site of the infinite system is ob- 
tained from the maximum eigenvalue A of the transfer 
matrix as F — — (T/2)lnA. Static quantities like spe- 
cific heat and susceptibilities are obtained from the free 
energy. The specific heat is calculated by the numeri- 
cal derivatives of F with respect to the temperature T. 
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The spin and charge susceptibilities are calculated by a 
shift of F under an applied magnetic field or chemical 
potential. 



G(iuj n ) or XAB(iw n ) by rational funcitons of frequency 
iu> n which are analytically continued to the real axis by 
iuj n — > lu + iS. 

The maximum entropy method is based on the spectral 
representations 



T(M=3) 



T(M=2) 





=2) 
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FIG. 5. Increasing of M of the quantum transfer natrix 
T. The Q n ,m is the indices of the new basis states. 

The calculation of the dynamic quantities requires ad- 
ditional steps. We first calculate correlation functions in 
the direction. This calculation requires good accuracy 
for the eigenvectors of the transfer matrix and \^ R ). 
Thus it is necessary to use the finite system algorithm of 
the DMRG. The Green's function in the (3 direction is 
obtained from the left and the right eigenvectors of the 
transfer matrix whose eigenvalue is the largest (See Fig. 

W- 

G{t 3 ) = -Tr{e-^c 4CT (r,) 4(0)}/Z 

= -(^1^(^)4(0)1^). (27) 

Similarly a local dynamic correlation function xab ) is 
obtained as 

X AB{Tj) = Tr{e~~ @ H Ai(rj) 5,(0)}/^ 

= (* l \Mtj) B t (0)\V R )- (28) 

By Fourier transformation, the Green's function and 
the the dynamic correlation function as functions of the 
imaginary frequencies are obtined as 



3 



XAB(iUn) 



(29) 
(30) 



G(r) = 
Xab(t) = 



1 + e-P" 



du>, 



1 e~™ 
-hn X AB(w) 1 _ e _p u> dw, 



(31) 
(32) 



with p(uj) = — — G(w + iS) being the density of state. 
Starting from a flat spectrum, this method finally finds 
optimal p(lo) and xab(^) that reproduce G(t) and 
Xab{t) best. 

The dynamical structure factor Sab{u) is related to 
the imaginary part of xab(^) through the fluctuation 
dissipation theorem, 
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FIG. 6. Schematic diagram of the calculation of the imag- 
inary time correlation function ^ h \ci a {rz) 4( T2 )l^ rH )' The 
Trotter number of this example is M = 3. 

In section III we use the standard zero-temperature 
DMRG method to study the ground-state properties of 
the paramagnetic metallic phase. In section IV after a 
brief discussion of the ground-state properties by using 
the zero-temperature DMRG, finite temperature prop- 
erties of the Kondo spin liquid phase will be discussed 
extensively by using the finite-T DMRG. 



where io n is the Matsubara frequency that is 7r(2n + 1)//? 
for fermionic operators and 27m/ ' (3 for bosonic operators. 

The real frequency Green's function and dynamic sus- 
ceptibility are obtained by the analytic continuation to 
the real frequency axis. We can use the Pade approxi- 
mations or the maximum entropy method for this pur- 
pose. The former method is based on the fittings of 



III. TOMONAGA-LUTTINGER LIQUID 
PROPERTIES OF THE PARAMAGNETIC 
METALLIC PHASE 

This section concerns with the paramagnetic phase 
away from half-filling, see Fig.|l|. Since there is no symme- 



G 



try breaking it is natural to consider that away from half- 
filling the translationally invariant Kondo lattice model 
is metallic. In one dimension it is well known that var- 
ious interacting metallic systems including the Hubbard 
model and the t-J model belong to the universality class 
of Tomonaga-Luttinger liquids jflXfl . Therefore, the first 
question concerning the paramagnetic metallic phase of 
the Kondo lattice model is whether it belongs to this class 
or not. 

The spin- 1/2 Tomonaga-Luttinger liquids have gapless 
charge and spin excitations. In one dimension the charge 
excitations are characterized by the velocity of the charge 
density v p and the correlation cxponcmt K p . Similarly, 
the spin excitations are characterized by the velocity of 
the spin density v a , but the correlation exponent in the 
spin sector is fixed by the SU(2) symmetry, K a = 1. 
The low-energy physics of a Tomonaga-Luttinger liquid 
is completely determined when these parameters are ob- 
tained. For example, the spin and charge susceptibilities 
are given by 



Xa = 



x P 



2K n 



nv 



(34) 
(35) 



Reflecting the gapless excitations, the density-density 
and spin-spin correlation functions show power-law de- 
cays where the exponents are determined by the correla- 
tion exponent, K p . The asymptotic forms of the density- 
density and spin-spin correlation functions are 

(n(x)n(0)) = K p /(ttx) 2 + Ax cos{2k F x)x- {1+K " ] 

+A 2 cos{Ak F x)x- iK » , (36) 
(S(x) • 5(0)) = 1/(ttx) 2 + B 1 cos(2k F x)x- (1+K >''> , (37) 

where k F — irp/2, with p being the density of charge 
carriers, is the Fermi momentum. The logarithmic cor- 
rections are neglected in Eqs.(|36|) and ( |37| ) |T^| . 

For the Hubbard model or the t-J model, the defini- 
tion of the density of carriers is straighforward. On the 
other hand, for the Kondo lattice model it is already 
questionable. When we take naively the conduction elec- 
trons as carriers, then the Fermi momentum is given by 
k F = k Fs = nn c /2. However, a different point of view is 
possible. Let us consider the Kondo lattice model as an 
effective Hamiltonian for the periodic Anderson model. 
For the latter, the density of carriers is the sum of the 
/ electon density and the conduction electron density. 
According to the Luttinger sum rule jl| the position of 
Fermi points do not change when the interaction is in- 
creased as long as the ground state remains paramag- 
netic. Therefore this property may be carried over to 
the Kondo lattice model and it would be also natural to 
assume k F = k F i = tt(1 + n c )/2 for the Kondo lattice 
model. 

Concerning the paramagnetic metallic phase, there are 
two basic questions: 



(1) Is it a Tomonaga-Luttinger liquid? 

(2) If it is the case, what is the size of the Fermi mo- 
mentum? Is it large, k F i = 7r(l + n c )/2 or small, 
k Fs = 7m c /2? The DMRG is a powerful method to ad- 
dress these questions. 

Let us define the ground-state energy in a given spin-5 
subspace for a finite system with L-sites by E g (L, N c , 5) 
where N c is the number of conduction electrons. In the 
following we will consider only even L and N c and thereby 
integer 5. The spin gap of a finite system is defined by 

A s (L)=E g (L,N c ,S = l)-E g (L,N c ,S = 0) . (38) 

Concerning the charge excitations we study the difference 
of the chemical potentials, p+ — /i_ which is given by 

2 M+ = E B {L, N c + 2,5 = 0) — E g (L, N c , 5 = 0), (39) 
2 M _ = E g {L, N c — 2,5 = 0) — E g (L, N c , 5 = 0). (40) 

Figure |?j shows (a) the spin excitation gap and (b) the 
difference of the chemical potentials as a function of in- 
verse of the system size. For the example the density of 
conduction electrons are fixed to n c = 2/3 fTsj] . Both 
quatities go to zero as 1/L — > 0, which means that the 
spin excitations and the charge excitations are gapless. 
Therefore it is most likely that the paramagnetic metallic 
phase of the Kondo lattice model belongs to the univer- 
sality class of the Tomonaga-Luttinger liquids. 

Now we will determine the parameters of the 
Tomonaga-Luttinger liquid. From the slope of the spin 
gap the velocity of the spin excitations are determined 
by 



A S (L) = v a ir/L, 



(41) 



since the lowest spin excitation of a finite system with 
the open boundary condition has the momentum ir/L. 
The difference of the chemical potentials is related with 
the charge susceptibility by 



M+ - M- = 



X P L 



(42) 



The values in Table I for v a and x P are determined from 
the slopes. The values of the spin susceptibility in the Ta- 
ble I is obtained from the spin velocity through Eq.(|34"|). 

To determine the charge velocity and the correlation 
exponent separately another independent measurement 
is necessary. The correlation exponent K p may be deter- 
mined from the density-density or the spin-spin correla- 
tion functions. However, determination of the exponent 
of a power law decay is the most difficult for any numeri- 
cal calculations. In order to determine K p we need to see 
the long-range behaviors of the correlation functions with 
sufficient accuracy. Instead of looking at the correlation 
functions we looked at the Friedel oscillations because 
the latter are numerically more reliable than the former 
P3|p^| . The reason for this is that the correlation func- 
tions are site off-diagonal, while the Friedel oscillations 
are site diagonal. 
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FIG. 7. (a) Size dependence of the difference of the chemi- 
cal potentials, /i+ — fi—, in the one dimensional Kondo lattice 
model, (b) Size dependence of the spin gap. L is the sys- 
tem size and the density of conduction electrons is fixed to 
n c — 2/3. The energy unit is t. Typical truncation errors in 
the DMRG calculations are 10 -4 . 



TABLE I. Luttinger liquid parameters of the one dimen- 
sional Kondo lattice model. The carrier density n c is 2/3. 
The energy unit is t. The errors are estimated from the am- 
biguity of the power law decay of the charge density Friedel 
oscillations. 







K P 






Vp 


Xp 


J/t = 





1 






1.73 


0.37 


J/t = 


1.5 


0.19 ± 0.03 






0.30 ± 0.06 


0.42 


J/t = 


1.8 


0.24 ± 0.02 


0.014 


46 


0.41 ± 0.06 


0.38 


J/t = 


2.0 


0.27 ± 0.02 


0.011 


56 


0.48 ± 0.06 


0.36 



The Friedel oscillations are density oscillations induced 
by a local perturbation. In a Tomonaga-Luttinger liquid, 
power law anomalies in correlation functions naturally 
reflect themselves in the Friedel oscillations. The usual 
Friedel oscillations induced by an impurity potential are 
given by 

Sp(x) ~ C\ cos(2k F x)x { ~ 1 ~ K "' )/2 + C 2 cos(4k F x)x~ 2K " 

(43) 



as a function of the distance x from the impurity [|15|- 
Analogously, spin density oscillations induced by a local 
magnetic field behave as 



cr(x) ~ D\ cos(2kpx)x 



(44) 



Thus, we can determine K p from the asymptotic form of 
the oscillations. It is worth to be noted that origin of 
the RKKY interaction may be traced back to the spin 
density oscillations induced by a localized spin. 

The charge density oscillations induced by the open 
boundary conditions are shown in Fig.|| for J = 1.5i and 
J = 2.5t at the density n c = 4/5 fll4| . The spin density 
oscillations induced by the local magnetic fields applied 
at the both ends with opposite directions are shown in 
Fig|] for the same set of parameters JlJ]. It is clearly 
seen that the dominant period of the charge density os- 
cillations is five sites, q — 2ir/5, while for the spin density 
oscillations ten sites, q = 7r/5. 

In the strong coupling limit of the Kondo lattice model, 
each conduction electron form a local singlet with the / 
spin on the same site. However away from half-filling 
these singlets can move in the lattice with the effective 
hopping matrix elements reduced by half. Equivalcntly, 
we can regard the unpaired / spins as mobile objects 
with the reduced hopping energy t/2 with its sign re- 
versed. Note that in the original model hopping matrix 
elements are defined by —t. Thus the effective model of 
the strong coupling limit is the t-model where the num- 
ber of carriers is L — N c and double occupancy of the 
carriers is prohibited. 

The charge response of the system is identical to the 
spinless Fermions where the Fermi point is given by 
7T— 7r(l— n c ) = im c . Therefore the induced charge density 
shows oscillations corresponding to 2im c which is equiv- 
alent to Akp s and 4fej?;. This analysis shows that in the 
strong coupling limit the amplitude of the Akp oscilla- 
tions dominates over the 2k p oscillations for the charge 
Friedel oscillations, Eq.(^). The period of five sites is 
naturally understood in this way and we have confirmed 
that for a weaker coupling the amplitude of the 2kp os- 
cillation develops. 
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FIG. 8. Charge density oscillations of the Kondo lattice 
model. The system size is 60 sites and the carrier density 
is n c — 4/5. The solid line and the broken line correspond 
to J = 2.5* and J — 1.5t, respectively. Typical truncation 
errors in the DMRG calculations are 1 x 10" 6 for J = 2.5* 
and 3 x 10" 6 for J = 1.5*. 
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FIG. 9. Spin density oscillations of the Kondo lattice 
model. The system size is 60 sites and the carrier density 
is n c = 4/5. The solid line and the broken line correspond 
to J = 2.5* and J = 1.5*, respectively. The strength of the 
local magnetic field h is 0.1*. Typical truncation errors in the 
DMRG calculations are 1 x 10~ 6 for J = 2.5* and 3 x 10~ 6 
for J = 1.5*. 

When we consider the spins of the unpaired / elec- 
trons, there remains a macroscopic 2 L ~ Nc -{o\d degener- 
acy. Concerning the spin sector the lifting of the de- 
generacy is essential. For the specific model where the 
degeneracy is lifted by the next nearest neighbor hop- 
pings it is shown analytically that the Fermi surface is 
big, kp — kpi |l8| . For a finite J the degeneracy of 
the Kondo lattice model is always lifted. The period of 



ten sites of the spin density Friedel oscillations shown in 
Fig. H indicate that the 2kp oscillations corresponding to 
2fcpi are actually observed. Figure [l(] show that the 2kFi 
oscillations corresponding to the large Fermi surface are 
always dominating in the paramagnetic phase for various 
coupling constants and various densities. 
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FIG. 10. Fourier components of the spin density Friedel 
oscillations. 

The Friedel oscillations obtained by the DMRG 
method clearly indicate that the Fermi surface of the 
Kondo lattice model is large. At the early stage, the 
conclusions of the bosonization studies were controver- 
sial. In the area of the paramagnetic metallic phase Fu- 
jimoto and Kawakami obtained the Tomonaga-Luttinger 
liquid with a large Fermi surface, while White and Affleck 
predicted a Luther-Emery liquid with a spin gap p9| , pfj[ . 
Later, it was argued that an additional direct Heisenberg 
coupling between the /^spins is necessary to stabilize the 
Luther-Emery liquid fl2lj . Furthermore, the existence of 
a gapless excitation with the momentum of 2kpi is shown 
rigorously by the Lieb-Schultz-Mattis construction . 

In order to obtain the correlation exponent K p , we 
used the slope of the envelope function of the charge den- 
sity oscillations, assuming that dominant component of 
the oscillations is the Akp oscillations. Figure |ll| shows 
K p thus determined for the exchange coupling constants 
from J = 4.0* to 1.5*. The density of conduction electons 
is fixed to n c — 2/3. K p is always smaller than 1/2 and 
monotonically decreases with decreasing J. The limiting 
value of K p — 1/2 in the strong coupling limit is easy to 
understand since the strong coupling limit of the Kondo 
lattice model is equivalent to the U — oo Hubbard model. 

The correlation exponent shows a small discontinuity 
at the boundary between the ferromagnetic and param- 
agnetic phases, J c = 2.4* for n c — 2/3. Below J c , the 
K p decreases faster and becomes lower than 1/3, which 
means that the long range behavior of the density-density 
correlation is governed by the Akp oscillations rather than 







the 2kp oscillations. With further decreasing J, the K p 
seems to cross the value 3 — 2y2 ~ 0.17. Since the ex- 
ponent of the power law anomaly of the momentum dis- 
tribution function is given by (K p + 1/K p — 2)/4, the 
power law anomaly is removed below this point and a 
clear Fermi surface can not be seen anymore. 

Through the study of the Friedel oscillations by the 
DMRG method, it has become clear that the paramag- 
netic metallic phase of the one dimensional Kondo lattice 
model is a Tomonaga-Luttinger liquid with a large Fermi 
surface. This Tomonaga-Luttinger liquid is unique in the 
sense that the K p is smaller than 1/2. The small K p may 
be attributed to the long range nature of effective inter- 
actions with strong retardation [ fL2|]23]| . Recently, obser- 
vation of the Friedel oscillations by the DMRG method 
is shown to be useful also for discussions of critical be- 
haviors of the Hubbard model IM. 
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/ 

FIG. 11. Correlation exponent K p estimated from the 
decay rate of the charge density Friedel oscillations. The er- 
rorbars are estimated from ambiguity of the power law fitting. 
n c = 2/3. J is in units of t. 



IV. KONDO SPIN LIQUID PHASE AT 
HALF-FILLING 



The half-filled KL model is always insulating in one- 
dimension. This conclusion was obtained by the exact 
diagonalization study [jj). Based on a finite size scaling 
it was shown that finite excitation gap remains for any 
finite J. The result has been confirmed by the DMRG 
method |25| l and later supported by the mapping to a 
non-linear sigma model and by the bosonization ap- 
proach [^7j . In this section we discuss basic properties of 
this insulating state. 

Concerning the insulating phase, the question we 
would like to address here is which characters distintin- 



guish the Kondo insulators from usual semiconductors. 
The most significant difference is that there are no gaps at 
high temperatures and they are induced as the tempera- 
ture is lowered. Furthermore the excitation gaps induced 
by the temperature are different depending on channels. 
The difference of the excitation gaps and more generally 
the difference in the temperature dependence of various 
excitation spectra are naturally reflected in temperature 
dependence of various thermodynamic quantities. 

Clearly the lowest excitation gap which is the spin gap 
for the Kondo insulator defines the smallest energy scale 
of the system. In ordinary band insulators the band gap 
defines the smallest energy scale which controlls not only 
the spin excitations but also the charge excitations. It is 
also interesting to compare the smallest energy scale of 
the Kondo lattice problem with that of the single impu- 
rity Kondo effect, namely the Kondo temperature Tk- It 
is well known that the low temperature properties of the 
impurity model is governed by the single energy scale of 

In the present section, first we will discuss the spin gap, 
the charge gap and the quasiparticle gap by the ground 
state DMRG. Then we will discuss temperature depen- 
dence of the spin susceptibility, the charge susceptibility 
and the specific heat by the finite-T DMRG. Temper- 
ature dependences of the single particle excitation spec- 
trum, the dynamic spin-spin correlation function and the 
charge-charge correlation function are also discussed by 
the finite-T DMRG, For the analytic continuation which 
is necessary to discuss the dynamic quantities, it is shown 
that the maximum entropy method is very useful [Eq-pOl . 



A. Spin, charge and quasiparticl gaps 

To understand the physics of the insulating state of 
the half-filled Kondo lattice model, it is instructive to 
consider the limit of strong exchange coupling J. In this 
limit every / spin together with a conduction electron 
form a local singlet on every site. To create spin excita- 
tions the minimum energy cost is J, which is the energy 
difference between the local spin singlet state and the 
local spin triplet state. On the other hand, creation of 
charge excitations requires the minimum energy of 3 J/2 
which corresponds to the energy cost for breaking two 
local singlets by transferring a conduction electron to a 
neighboring site. 

The excitation gaps monotonically decrease with de- 
creasing exchange constant, but they do not vanish at 
any finite value of J. Particularly, the weak coupling 
limit J -C t is interesting. In this regime the KL model 
is equivalent to the periodic Anderson model with strong 
Coulomb repulsion in the / orbitals. The salient fea- 
ture of the strong Coulomb interaction in the periodic 
Anderson model appears in the diverging ratio between 
the charge and spin gaps. The limit of J — is singu- 
lar where the conduction electrons and the / spins are 
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decoupled, and both the spin and charge gaps vanish. 

For the discussion of the gaps, we take into account also 
the Coulomb interaction between the conduction elec- 
trons. Since the spin and charge gaps are tiny in the weak 
coupling regime, it is no more justified to neglect the 
Coulomb interaction between the conduction electrons. 
This Coulomb interaction suppresses the double occupa- 
tion of conduction electrons, which eventually leads to 
the formation of local magnetic moments of the conduc- 
tion electrons. Therefore the effect of the Coulomb in- 
teraction on the spin and charge gaps of the KL model 
sheds lights on the nature of the gap formation in Kondo 
insulators. 

The model we consider in this subsection is the follow- 
ing one-dimensional KL model with Coulomb interaction 
between the conduction electrons U c : 

H = -*E( c *Wi. + H - c -) + JEW 

ia ij-L 

+ u c £>W - - \)- ( 45 ) 

i 

The Coulomb interaction is represented in the last term. 
In this section we consider the case of half-filling where 
the total number of conduction electrons is equal to the 
number of lattice sites L: N c = J2ia c \a c i<? = L. This 
Hamiltonian is reduced to the Hubbard model in the limit 
of J — > 0, and to the usual KL model for U c = 0. 

In the impurity Kondo model all low temperature 
properties are scaled by the single energy scale Tr- ~ 
.Dexp (—jj), where p is the density of states of the con- 
duction band at the Fermi level and is given by ^ in 
one dimension. In contrast to the single impurity Kondo 
model, the KL model has many / spins which are cou- 
pled through the conduction electrons. A basic question 
of the lattice problem is how the intersite correlations ap- 
pear in the energy scale. The simplest extension of the 
form of Tk may be an inclusion of an enhancement factor 
in the exponent: The spin gap is expected to behave as 



A S (L) = A s (oo) + PL- 2 + 0{L- 4 ). 



(47) 



A s oc exp ( ) 

apJ 



(46) 



where a is the enhancement factor. The Gutzwiller ap- 
proximation predicts an enhancement factor a — 2 [5IJ. 
Tsunetsugu el al. have estimated that the enhancement 
factor in one dimension is in the range of 1 < a < 5/4 
by using a finite size scaling for the results obtained by 
the exact diagonalization j{| . In the following we present 
the results on the enhancement factor obtainedy by the 
DMRG. 

The spin gap is obtained from the difference of the 
ground-state energies in the subspaces of total S z being 
zero and one, Eq.(H|); the SU(2) symmetry in the spin 
space guarantees the energy difference is the same as the 
spin gap in the subspace of zero total S z . The spin gap 
of the bulk system is estimated from the following scaling 
function: 



The obtained spin gaps are plotted in Fig. [12| in log- 
arithmic scale as a function of 1/J. The results are 
obtained by the extrapolation to the bulk limit using 
data of L = 6, 8, 12, 18, 24, 40. The DMRG calculations 
were done by using the finite system algorithm with open 
boundary conditions keeping up to 300 states for each 
block. The enhancement factor is obtained from the slope 
in the figure, and determined to be a = 1.4(1) for U c = 0. 
There are some uncertainties in the extrapolation to the 
bulk limit for tiny gaps. However within the present accu- 
racy we do not observe any indication of the logarithmic 
correction to the exponent which was predicted by the 
semiclassical approach 26 . 
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FIG. 12. Spin gap of the half-filled one-dimensional Kondo 
lattice model with Coulomb interaction. The thick curve rep- 
resents the result of the perturbation theory in terms of t/J 
for U c = lOt. Typical truncation error in the DMRG calcu- 
lation is 10 -6 for J = 1. Errorbars are estimated from L _1 
and L~ 2 scalings. Gap energies, exchange constant J, and 
Coulomb interaction U c axe in units of t. 

Now we consider the effect of the Coulomb interac- 
tion. In the weak coupling region it is natural to extend 
the form of Eq. (46) to finite U c allowing J7 c -dependence 
of the exponent. Indeed the numerical data are nicely 
fitted by this form as shown in Fig. [l2|. The obtained 
?7 c -dependence is shown in Fig. which indicates that 
ct(U c ) increases with increasing U c and the asymptotic 
behavior is linear in U c . The KL model with the Coulomb 
interaction is mapped to a Heisenberg chain coupled with 
the localized / spins in the limit of U c /t — > oo. The lin- 
ear J7 c -dependence of the exponent a = 0.78U c /t + 0.7 
in Fig. [l3| means that the spin gap of the effective spin 
model behaves as A s ~ exp (— 2 J Q s/ J) with J c g = 4t 2 /U c 
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beeing the effective coupling of the Heisenberg chain. In 
order to check this form we have analyzed the numeri- 
cal data for the spin system obtained by Igarashi et al. 
p2| , and found a good coincidence. Thus we conclude 
that the enhancement factor a increases monotonically 
with increasing U c . This is natural since the origin of 
the spin gap is the singlet binding between the localized 
spins and the conduction electons for which the Coulomb 
interaction helps by suppressing the double occupancy. 

In contrast to the single impurity Kondo model, the 
KL model has the second energy scale that characterizes 
the charge excitations at low temperatures. The charge 
excitations keep spin quantum numbers, and the charge 
gap is defined by the difference of the lowest energy in 
the subspace of N c = L and N c = L + 2: E g (L,N c — 
L + 2,5 = 0) - E g (L, N c = L,S = 0). Owing to the 
hidden SU(2) symmetry in the charge space, the energy 
difference is the same as the charge excitation gap in the 
subspace of the fixed number of electrons N c — L [R3| . 




Uc 

FIG. 13. [/c-dependence of the exponent of the spin gap. 
Coulomb interaction U c is in units of t. 
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shows the charge gap obtained by the ex- 
trapolation to the infinite system. The results for J = 
are known as the Hubbared gap of the one-dimensional 
Hubared model which is exactly solved by the Bethe 
Ansatz |$4|. The asymptotic forms of the charge gap 
are given by A c oc \fU c t exp (—^7^) for small U c , and by 
A c ex U c — 4t for large U c . The results obtained for finite 
J are consistent with the exact ones which are denoted 
by the crosses on the vertical axis. 

For U c = the charge gap is linear in J in the small 
J/t limit. As is shown by the exact diagonalization study 
the charge gap is much bigger than the spin gap in the 
weak coupling regime |]33f . It implies that the correlation 
length for the spin degrees of freedom is much longer 



than the charge correlation length. Therefore for the 
discussion of the charge gap it is justified to assume that 
the spin-spin correlation length is infinitely long. Under 
the assumption of the infinite spin correlation length, the 
charge gap is calculated as 



J 

A c = -. 
2 



(48) 



From Fig. [Ij we also find that the charge gap increases 
with increasing Coulomb interaction U c . 
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FIG. 14. Charge gap of the half-filled one-dimensional 
Kondo lattice model with Coulomb interaction. Results on 
the vertical axis are obtained from the exact solution by 
Lieb-Wu. Typical truncation errors in the DMRG calcula- 
tion are 10 -6 for J = 1 and 10~ 4 for J = 0.2, which are 
dominant source of numerical errors since the finite size scal- 
ing, Eq. ([47]), is well obeyed. Gap energies, exchange constant 
J, and Coulomb interaction U c are in units of t. 

The charge excitations are created by adding two ad- 
ditional electrons keeping spin quantum numbers fixed. 
When we put single electron in the ground state, then 
a quasiparticle excitation is made. Here we consider the 
relation between the charge gap A c and the quasiparti- 
cle gap A qp which is defined by E g (L, N c = L ± 1, 5 = 
±1/2) - E g (L, N c = L,S = 0). In the strong coupling 
limit, J/t — > 00, it is evident that the charge gap is twice 
the quasiparticle gap owing to the SU(2) symmetry in 
the charge space. In the second order perturbation in 
t/J, one can show that the interaction between the two 
additional electrons is repulsive, leading to only a phase 
shift. Therefore the charge gap in the bulk limit is twice 
the quasiparticle gap A qp : 



2A„ 



(49) 



A similar argument is also valid for the periodic Ander- 
son model I33J . Validity of this relation is checked by 
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the DMRG calculation in the entire region of the ex- 
change constant J. Concerning the spin gap, the lowest 
spin excitation may be considered as a bound state of a 
quasielectron and a quasiholc. 



B. Susceptibilities at finite temperatures 

The spin and charge gaps determined at zero tempera- 
ture are very different in the weak coupling regime. The 
spin gap is exponentially small, while the charge gap is 
proportional to J. The large charge gap originates from 
the staggering internal magnetic fields induced by the 
long correlation length of the / spins. At finite tem- 
peratures, however, the spin correlations are subject to 
the thermal fluctuations. When temperature becomes 
comparable to the spin gap, the spin correlation length 
gets smaller and the whole electronic states including the 
charge excitation spectrum are reconstructed. In this 
section we study such an interplay between the spin and 
charge excitations at finite temperatures by looking at 
the thermodynamic quantities. In what follows we con- 
sider the original KL model neglecting the Coulomb in- 
teraction between the conduction electrons. 

In order to calculate thermodynamic quantities we use 
the fmite-T DMRG discussed in Sec. II @||. In this 
method the free energy is obtained from the maximum 
eigenvalue of the quantum transfer matrix. The spin and 
charge susceptibilities arc obtained from the derivatives 
of the free energy with respect to external magnetic field 
or chemical potential. The calculations are performed 
by the infinite system algorithm keeping 40 states per 
block. The truncation errors in the calculation are typ- 
ically 10 -3 and at the lowest temperature 10 -2 for the 
Trotter number M = 50. 
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FIG. 15. Spin susceptibility of the half-filled 

one- dimensional Kondo lattice model. The truncation errors 
in the finite- T DMRG calculation are typically 10 -3 and at 
the lowest temperature 10 -2 . 



uniform spin susceptibility. The spin susceptibility is ob- 
tained from the change of the free energy by a small 
magnetic field h: 6F — Xsh 2 /2 The results for J/t — 
0, 1.0, 1.2, 1.6 and 2.4 are shown in Fig. ||. 

When J/t = 0, the localized spins and conduction elec- 
trons are uncorrelated. The susceptibility is given by the 
sum of the Curie term due to the free / spins and the 
Pauli term of the free conduction electrons. The contri- 
bution of the Pauli susceptibility of the conduction elec- 
trons shown by the dashed line in Fig. [l5| is relatively 
small, and the total susceptibility for J/t — is domi- 
nated by the Curie term. 

For finite J, the low temperature part of \ s sharply 
drops with decreasing the temperature. This drastic 
change is due to the appearance of a low energy scale 
for the spin sector. The spin gap of 0.08£ for J/t = 1.0 is 
consistent with the characteristic temperature at which 
Xs starts to decrease deviating from the Curie low. 

In order to determine the energy scale at low tempera- 
tures we estimate the activation energy by fitting the sus- 
ceptibility with an exponential form. The estimated ac- 
tivation energy for the spin susceptibility is summarized 
in Table || for J/t = 1.6, 2.4 and 3.0. Compared with the 
quasiparticle gap and the spin gap, both of which are re- 
sponsible for the magnetic excitations, we conclude that 
the lower one of them determines the low-temperature 
energy scale of the spin susceptibility. This is consistent 
with the general form of susceptibility which is written 
as 



m 

Z = E< 



-PE„ 



(50) 
(51) 



The point is that Eqs.(50) and ( |5l| ) apply for both the 
canonical and the grand canonical ensembles by properly 
defining the states |m). In the thermodynamic limit the 
susceptibilities by the two ensembles should give the same 
answer. From this consideration it is concluded that the 
smaller one of the spin gap and the quasiparticle gap 
determines the low temperature energy scale. For the 
case of small exchange coupling ( J/t <C 1) the spin gap 
is smaller than the quasi-particle gap and thus the low 
temperature energy scale of Xs is determined by the spin 
gap. 

In order to see the effect of thermal fluctuations of 
/ spins to the charge excitations, we next calculate the 
charge susceptibility Xc- The Xc is obtained from the 
change of the free energy due to a small shift of chemical 
potential [i, SF = x c /i 2 /2. In the present calculation we 
use the fact that the chemical potential is zero at half- 
filling owing to the SO(4) symmetry of the model. |33j ] 
The results for J/t — 0, 1.0, 1.2, 1.6 and 2.4 are shown in 
Fig. M 



We first consider the temperature dependence of the 
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FIG. 16. Charge susceptibility of the half-filled 

one- dimensional Kondo lattice model. 

For J/t = 0, Xc does not show diverging behavior at 
low temperatures in contrast to Xs- In the limit of T = 0, 
Xc is equal to the density of state of conduction electrons 
which is 1/irt including the two spin directions. This is 
expected since the charge degrees of freedom is governed 
by the conduction electrons. Since there is no interac- 
tion, Xc/4 is equal to the spin susceptibility of the free 
conduction electrons. The slight increase in Xc in the 
low-temperature region is a characteristic feature of the 
one-dimensional system where the density of states di- 
verges at the band edges. 

A finite value of J produces a sharp drop in Xc at low 
temperatures. Similarly to Xs, this drop is due to the 
appearance of low energy scale A c for the charge sector. 
The energy scale is determined by the activation energy 
for the charge susceptibility. By fitting Xc with an ex- 
ponential form, the activation energy A Xc is obtained as 
listed in Table [0] for J/t = 1.6, 2.4 and 3.6. From this ta- 
ble it is concluded that the quasiparticle gap determines 
the low-temperature energy scale of the charge suscepti- 
bility. 

TABLE II. Activation energy obtained from the spin and 
charge susceptibilities, A Xs and A Xc , of the one- dimensional 
Kondo lattice model at half-filling. The quasiparticle gap A qp 
and the spin gap A s are obtained by the zero-temperature 
DMRG. The charge gap is twice the quasiparticle gap; 
Ac — 2Aqtj 







Ax, A 


Axe A 


A a /t 


A qp /t 


J/t = 


1.0 






0.08 


0.36 


J/t = 


1.2 






0.16 


0.47 


J/t = 


1.6 


0.45 ±0.1 


0.6 ±0.1 


0.4 


0.7 


J/t = 


2.4 


1.2 ±0.1 


1.0 ±0.1 


1.1 


1.1 


J/t = 


3.0 


1.6 ±0.1 


1.4 ±0.1 


1.8 


1.5 



Although the qusiparticle gap determine the exponen- 
tial temperature dependence at low temperatures, it is 
not the single energy scale for Xc- It may be best under- 
stood by looking at Xc for J/t = 1.0. A sharp decrease 
of Xc is seen at around T ~ O.lt which is much smaller 
than the quasiparticle gap 0.36i but rather close to the 
value of the spin gap 0.08i. This fact suggests that the 
whole electronic states are reconstructed when temper- 
ature is raised up to the spin gap. We will discuss this 
aspect in more detail in connection with the temperature 
dependence of various excitation spectra. 



C. Specific heat 

In order to see how the entropy of the system is re- 
leased, we next calculate specific heat. The specific heat 
is calculated from the second derivative of the free energy; 
C = -Td 2 F/dT 2 . The results for J/t = 0, 1.0, 1.2, 1.6 
and 2.4 are shown in Fig. [l?]. 

At J/t = the specific heat of this model is given by 
the sum of the delta function at T = that originates 
from the free localized spins and the specific heat of free 
conduction electrons. For finite J they are combined to 
form a two-peak structure. The peak at higher temper- 
atures is almost independent of the exchange constant 
and similar to the specific heat of free conduction elec- 
trons. Thus the structure at higher temperatures may 
be understood as the band structure effect of the one di- 
mensional conduction electrons. In contrast to the higher 
temperature structure, the structure at lower tempera- 
tures strongly changes its form with J. The peak shifts 
towards higher temperatures and becomes broader with 
increasing J. 

With further increasing the exchange coupling, the 
spin gap becomes comparable to the hopping matrix el- 
ement t. In this situation various energy scales are not 
distinguishable and the specific heat possesses a single 
peak structure as shown for J/t = 2.4. 
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FIG. 17. Specific heat of the half-filled one-dimensional 
Kondo lattice model. 



D. Dynamic properties 

Dynamic properties of the half-filled KL model show 
clear features characteristic to strongly correlated insula- 
tors. Unusual temperature dependence of the excitation 
spectra is one of the most important features of the in- 
teracting systems. Actually the behaviors of the static 
susceptibilities discussed in the preceding subsection in- 
dicates that the excitation gaps develop at low tempera- 
tures. 

In this section we calculate the dynamic spin and 
charge structure factors, S(w) and N(ui), and the density 
of states, p(to). By looking at their temperature depen- 
dence we can study temperature evolution of the excita- 
tion gaps and the relations among dynamic quantities. 
To obtain the dynamic quantities we first calculate the 
correlation functions in the imaginary time direction. As 
we have discussed in section II, the correlation functions 
in the imaginary time direction are directly calculated 
from the left and right eigenvectors obtained by applying 
the finite- T DMRG to the transfer matrix. 

Examples of the single particle Green's function as a 
function of the imaginary time calculated by the finite-T 
DMRG are shown in Fig. [t8]. To obtain the spectral 
functions we first Fourier transform the imaginary time 
correlation functions and then need to perform analytic 
continuation from the imaginary frequency axis to the 
real frequency axis. One straightforward method for the 
analytic continuation is to use the Pade approximations. 
Since the DMRG calculation yields no statistical errors 
the Pade approximations show good convergence in many 
cases. However, it is still difficult to obtain the spectral 
functions by the Pade approximations in a stable manner 
when a spectrum has a nearly singular form. The reason 
is that the Pade approximations use rational functions of 
Matsubara frequencies iuj n . 

Even when a spectrum has a nearly singular form, the 
maximum entropy method still works well. An advantage 
of this method is that we can explicitly use the symme- 
tries and the positiveness of the spectral function. The 
results obtained by the two method are compared in Fig. 
[l9| . At high temperatures two results coincide with each 
other, but with decreasing the temperature the conver- 
gence of the Pade approximations becomes worse due to 
the growing singularity in the spectral function in the low 
frequency region. The results by the maximum entropy 
method are stable even at low temperatures. Thus in 
the following we employ the maximum entropy method 
to study temperature evolution of the dynamic correla- 
tion functions iBal. 
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FIG. 18. Imaginary time correlation functions, G(r), 
of the half-filled one-dimensional Kondo lattice model: 
J/t — 1.6. The Trotter number M = 60. Inset shows the 
results for diffenent number of state kept m in the DMRG 
calculation. The truncation errors in the DMRG calculation 
are 2 x 10" 3 for m = 26 and 3 x 10~ 4 for m = 40. 
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FIG. 19. Quasiparticle density of states, p(u), obtained 
from the G(r) by the Pade approximations and the maximum 
entropy method (MEM). The Trotter number M = 60. 

The quasiparticle density of states obtained for J/t = 
1.6 at temperatures T/t = 0.1, 0.14, 0.2, 0.25, 0.3, 0.6 is 
shown in Fig. |2^. Existence of the quasiparticle gap is 
seen as a clear dip structure around u = 0. At low tem- 
peratures sharp peaks appear at u) = ±A gp separated 
from the higher frequency part of the spectral weight. 
Sharpness of the peaks suggests the formation of the 
heavy quasiparticle bands at the gap edges. The high 
frequency part of the spectral weight extends to the re- 
gion far from the edge of the free conduction band u> = 2i, 
which shows significance of the multiple excitations ac- 
companied with the quasiparticle excitations. 
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FIG. 20. Quasiparticle density of states, p(uj), of the 
half-filled one-dimensional Kondo lattice model: J/t — 1.6. 
The Trotter number M — 60, and the number of states kept 
m = 40. 
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FIG. 21. Dynamic spin structure factor of the / spins, 
Sf(ui), of the half-filled one-dimensional Kondo lattice model: 
J/t = 1.6. 



As the temperature is increased, the peak at the 
threshold gets broadened and at temperatures T ~ A s 
the peak and the accompanied dip between the low and 
high frequency parts completely disappear. This result 
shows that although the sharp peaks of the quasiparticle 
density of states are located at the frequency u> — ±A 9P , 
it disappears around the temperature T ~ A s which is 
much lower than T ~ A gp . 

In order to see what happens at the temperatures T ~ 
A s , we next consider the / spin dynamic structure factor 
Sf(u). The calculated Sf(w) are presented in Fig. 53. 
At the lowest temperature the spin gap is clearly seen 
with a sharp peak at the gap edge. This characteristic 
peak has the most of the spectral weight, which shows 
concentrated / spin excitations at the energy scale of 
A s . There is a broad peak in the higher frequency side. 
As is shown later, a similar structure and temperature 
dependence appear in the dynamic spin structure factor 
of the conduction electrons S c (uj). Thus we conclude that 
through the exchange coupling excitations of the / spins 
are mixed with those of conduction spins, which yields 
this broad peak in higher frequency part of Sf(uj). 

With increasing the temperature, the peak structure 
at u> = A s becomes broad and the spectral intensity 
increases around the zero frequency \lu\ < A s . At the 
temperatures T ~ A s , the peak position of the spectrum 
shifts to the zero frequency, and the peak height becomes 
almost temperature independent. The spectral intensity 
at the zero frequency is directly related to the NMR re- 
laxation rate \/T\. Hence the present results show that 
1 /T\ is nearly temperature independent at high temper- 
atures and drastically decreases with decreasing temper- 
ature below the characteristic temperature of the order 
A.,. 
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FIG. 22. Dynamic spin structure factor of the conduction 
electrons, S c (ui), of the half- filled one-dimensional Kondo lat- 
tice model: J/t = 1.6. 

The dynamic spin structure factor for the conduction 
electrons, S c (<^), is shown in Fig. At low temper- 
atures, S c (u>) has two peaks. The peak in the low fre- 
quency side is located at the energy of A s similarly to 
Sf(u>). This peak corresponds to the spin excitations 
of the singlet bound states of conduction electrons with 
/ spins to triplet states. The high frequency peak is 
located slightly above the charge gap, that corresponds 
to the spin excitations due to quasiparticles. With in- 
creasing the temperature both peaks lose their intensity, 
and above the temperature T ~ A s the low frequency 
peak structure disappears. The spectrum at the high 
frequency side becomes similar to that of the dynamic 
charge structure factor N(ui). This means that high fre- 
quency excitations are dominated by the quasiparticle 



1G 



excitations of almost free conduction electrons and thus 
the relation S c (uj) = N c (uj)/A is satisfied approximately. 



0.15 




FIG. 23. Dynamic charge structure factor of conduction 
electrons, N c (ui), of the half-filled one-dimensional Kondo lat- 
tice model: J/t = 1.6. 

The dynamic charge structure factor N(uj) is shown 
in Fig. B3l At the lowest temperature two clear peaks 
appear, a smaller peak at u> ~ and a bigger one at A c . 
These two peaks originate from the sharp peak struc- 
ture in p(oj) at uj = ±A gp . The excitations of thermally 
populated quasiparticles within the sharp peak in p(oj) 
contribute to the peak at uj = 0, while the excitations 
between the peaks in p(u>) make the peak in N(uj) at 
uj ~ A c . With increasing temperature, increased num- 
ber of the thermally populated quasiparticles enhances 
the peak at uj = 0, but at the temperature T ~ A s the 
peak structure is completely smeared out, which reflects 
the disappearance of the peak in p(uj). The gap structure 
of N(uj) and the energy scale of A c become unclear at 
temperatures much smaller than A c . 

Dynamic quatities studied by the finite-T DMRG have 
revealed the many-body nature of the gap formation in 
the Kondo insulators. The difference among excitation 
gaps depending on channels is a characteristic feature of 
the Kondo insulators compared with the ordinary band 
insulators. The temperature induced gap formation for 
the single particle density of states is another clear ev- 
idence of the many-body feature. Even at a fixed tem- 
perature a renormalizcd band picture fails to capture the 
essential physics of the strongly correlated insulators. A 
typical example is that the two-body excitation spectrum 
N(lu) is very different from a convolution of the one-body 
excitation spectrum p(uj). In the Kondo insulators there 
are several low energy scales, corresponding to the spin 
gap, the quasiparticle gap and the charge gap. Among 
them the lowest one, the spin gap, plays a special role. 
At higher temperatures than the spin gap, the excitation 
spectra in the charge sector are also modified strongly. It 
means that the whole electronic states are reconstructed 



above the tempeature corresponding to the lowest energy 
scale. 



V. SUMMARY AND DISCUSSIONS 

In this review we have discussed Tomonaga-Luttinger 
liquid properties of the one-dimensional Kondo lattice 
model away from half-filling. In particular, the large 
Fermi surface is concluded in the ground state by in- 
vestigating the spin and charge Friedel oscillations. At 
half-filling of the one-dimensional Kondo lattice model 
the ground state is always an incompressible spin liq- 
uid phase. Studies on the dynamic correlation functions 
have revealed many-body nature of this insulating phase 
in several ways. 

These developments have been achieved by applying 
the density matrix renormalization group method either 
to the Hamiltonian itself or to the quantum transfer ma- 
trix. In the problem of Kondo lattice model there ap- 
pear small energy scales at low temperatures. It implies 
that the correlation lengths for various quantities are rel- 
atively long and therefore we need sufficiently long sys- 
tems to observe intrinsic properties. On the other hand, 
there are 8 states per site in the Kondo lattice problem. 
Of course for a quantum spin-1/2 chain there are only 
two states per site. Exact diagonalization studies can do 
a good job for the latter but only poor one for the former. 
In this situation the DMRG shows its full advantage for 
the Kondo lattice model. 

We would like to stress that now we can calculate dy- 
namic quantities at finite temperatures by applying the 
finite-T DMRG to the quantum transfer matrix. This 
method is free from statistical errors and the truncation 
errors are the only source of numerical errors. Therefore 
much better accuracy is obtained for the imaginary time 
data, from which corresponding spectral function may be 
obtained reliably through the maximum entropy method. 
Another advantage of the finite-T DMRG compared with 
the quantum Monte Carlo simulations is that we do not 
have the negative sign problem for any quantum systems. 

Generally speaking, more elaborate calculations are re- 
quired for the finite-T DMRG compared with the zero- 
temperature DMRG. Finite temperature properties, both 
static and dynamic, of the Tomonaga-Luttinger liquid 
phase are of great interest and studies in this direction 
are now being in progress. In near future it will become 
possible to address these questions. Investigation of a 
completely Fermionic model for heavy Fermions, for ex- 
ample the periodic Anderson model, is also left for future. 
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